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Abstract 

We discuss systematic extensions of the standard (Stormer-Verlet) splitting method 
for differential equations of Hamiltonian mechanics, with relative accuracy of or- 
der T 2 for a timestep of length x, to higher orders in x. We present some splitting 
schemes, with all intermediate timesteps real and positive, which increase the rel- 
ative accuracy to order x N (for N = A, 6, and 8) for a large class of Hamiltonian 
systems. 
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1. Introduction 

The Hamilton equations of motion constitute a system of ordinary first order 
differential equations, 

■a dH . dH 

q df' Pa = ~d^' a=l,...,N, (1) 

where the ' denotes differentiation with respect to time t, and H = H(q,p). They 
can be viewed as the characteristic equations of the partial differential equation 

-^p(q,p;t)=Sfp(q,p;t), (2) 



with the first order differential operator, 

y dH d 

a =i d Pa dq a dq a dp, 
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generating a flow on phase space. If H does not depend explicitly on t, a formal 
solution of ([2]) is 

p(q,p;t) = e t ^p(q,p;0). (4) 

In most cases this expression remain just formal, but one may often split the 
Hamiltonian into two parts, H = H{+ H2, with a corresponding splitting = 
J?fi + J?2 such that the flows generated by ££\ and J??2 are separately integrable. 
One may then use the Cambell-Baker-Hausdorff formula to approximate the flow 
generated by Jz? . By doing this in a symmetric way one obtains the Strang splitting 
formula EM 

exp exp[T^i]exp[i^ 2 ] 

= exp [TJSf + ir 3 [2JSfj + J2f 2 , [^1 ,-%]] + -•-] , (5) 

which shows that time stepping this expression with a timestep % provides an ap- 
proximation with relative accuracy of order t , exactly preserving the symplectic 
property of the flow. 

This corresponds to the symplectic splitting scheme of iterating the process of 
solving 

for a timestep \-%, 





dH 2 




dH 2 


q a = 


dp a ' 


Pa = " 


dq a 




dHi 




dHi 


q a = 


dp a ' 


Pa = - 


dq a 




dH 2 




dH 2 


q a = 


dp a ' 


Pa = " 


dq a 



for a timestep T, (6) 
1 

for a timestep -T. 

Here the last part of one iteration may be combined with the first part of the next, 
unless one deals with time dependent systems or wants to register the state of the 
system at the intermediate times. 

There have been several approaches to construct integration schemes which 
are of higher order in T, while maintaining the exact symplectic nature of the evo- 
lution. Accessible reviews of such approaches have f.i. been given by Yoshida 0, 
McLachan et. al. |0, and Blanes et.al. 0. Neri [6] has provided the general idea 
to construct symplectic integrators for Hamiltonian systems. Forest and Ruth 
discussed an explicit fourth order method for the integration of Hamiltonian equa- 
tions for the simplest non-trivial case. Suzuki [8J presented the idea of how recur- 
sive construction of successive approximants may be extended to other methods. 
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Many of the higher order symplectic splitting methods involve an extention of 
equation ([5]) to an expression of the form 



fj exp [ax^ 2 ] exp [ekcSfi] = exp + Sf 2 ) + 



(7) 



as discussed by Yoshida I0. It was noted that if one uses a symmetric integrator, 
such that 

5(t) = exp [t(^i +^ 2 ) + T*R t + 0(T e+2 ) 
for some generator R(, then 

S(jcit) S{x q x) S{xi%) = exp Uxq + 2x 1 )t(£ ? 1 + ^ 2 ) + (jcJ+2*{ + ff(r e+2 ) 



Hence, by choosing 



(8) 
(9) 



xq + 2x{ = 1, 
x f + 2x[ = 0, 

one increases the order of the scheme by two or more. However, equations 
have real solutions only if either xq or x\ is negative. In fact, it has been proven 
(cf. Sheng 0H, Suzuki 0J], Goldman and Kaper [12], Blanes and Casas 031) that 
all schemes of the form <[7j) require at least one c, < 0, and at least one d{ < 0. For 
equations invariant under time reversal, which is often the case for Hamiltonian 
systems, this may not be a big obstacle (although it seems like an inefficient way 
of integrating equations forward in time). 

Worse, if one wants to use the same code to solve parabolic equations (like 
a heat type equation, —d t u(t,x) = [— A + V(x)] u(t,x), instead of a Schrodinger 
equation, id t y/(t,x) = [—A + V(x)] y/(t,x), (which formally corresponds to re- 
placing t by —it in the Schrodinger equation) negative timesteps may have a dis- 
astrous effect on numerical stability due to exponentially growing errors. Castella 
et. al. |[T4ll have proposed to use complex solutions of equations ([8j|9]). It is possi- 
ble to find solutions where all timesteps have a positive real part. This can stabilize 
the scheme, but at the cost of working with complex variables. 

In this paper we investigate a different approach, based on our [15J observation 
that the operators Jz?i, J£ 2 of each step of a splitting scheme don't need to be 
exactly the same as those in the sum J£ = J£\ + Jtf 2 . Instead, our approach is to 
construct T-dependent operators J£\ , ££ 2 such that 



exp [jT^ 2 ] exp [tJS?J exp [\x^ 2 ] = exp \%{%\ +^ 2 ) + G{x 



(10) 
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For Hamiltonians of the form, 

H(q,p) = ^p T Mp + V(q), (11) 

with M a symmetric positive definite matrix (the inverse mass matrix), we have 
constructed explicit expansions, 

(N-2)/2 

X = ^+ £ x 2k ^ k \ (12) 

k=\ 

for N = 4, 6, 8, leading to schemes with global error of order x N . We denote N the 
order of these schemes. Since the operators Jz^ generate flows exp [tJz^-] which 
are modifications of those generated by J^, we refer to such flows as modified 
integrators. Chartier et. al. lfT6ll have labeled such schemes as modified differential 
equations. Thus, the ATth order scheme is constructed to generate the same flow 
as 

JSf = JSf - [^ N/2) + 4 W2) ) + ^(t w+2 ), (13) 

when averaged over timesteps. I.e., we use modified integrators to generate the 
unmodified flow better. 

One possible restriction on the class of available splitting schemes is the re- 
quirement that both of the flows exp [jTJ?^] and exp should be explicitly 
computable. We have relaxed this requirement by demanding both flows to be 
efficiently computable: I.e., each (short) timestep must be possible to integrate nu- 
merically sufficiently fast, while preserving the symplectic structure to sufficient 
numerical precision. 

We have implemented this through the systematic construction of a generating 
function £f for an exact symplecting transformation which reproduces the flow 
over each finite timestep x to sufficient accuracy. For nonlinear dynamics this 
results in implicit formulas for the evolution 

: {q(t),p(t)} -+ {q(t + *),p(t + T)} ■ (14) 

This leads to schemes where each kick step is explicitly computable, while each 
move step actually involve a small implicitly defined push followed by an explic- 
itly computable move. Here kick and push involve a change in momenta p with 
fixed positions q; a move involves a change in positions q with fixed momenta p. 

The rest of this paper is organized as follows: In section [2] we demonstrate 
the basic idea of the proposed methods on linear systems. Next we develop the 
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general theory, valid for separable Hamiltonians ( fTT| ), in section |3j Here we con- 
struct the operators Jzf; explicitly; more precisely the Hamiltonians T21 and V^, 
cf. equation (30) corresponding to these operators. We focus our discussion on the 
numerical implementation of these methods in section [5} together with investiga- 
tions of how they work in practice. We have tested these methods on anharmonic 
oscillators and Fermi-Pasta-Ulam-Tsingou type problems (named as suggested by 
Dauxois ifTTlO . We close this paper with some brief concluding remarks in sec- 
tion |6] 



2. Linear Systems 

2.1. Single Harmonic Oscillator 

For a simple illustration of our idea consider the Hamiltonian 

H(q lP )= l -(p 2 + q 2 ), (15) 



whose exact evolution over a time interval x is 

<? e \ / cost sinr 
n e I V — sin T cos X 



(16) 



Compare this with the process of first evolving the system with the Hamiltonian 
//kick = \kq 2 for a time \x, followed by an evolution with //move = \mp 2 for a 
time x, and ending with an evolution by //kick = \kq 2 for a time \x. One such 
combination gives 



/s i j 1 — jmkz 2 mx 



^ 1 ' -(\ — \kmx 2 )kx \ — \kmx 




(17) 



We note that by choosing 



sinr 1 1 2 1 4 1 6 

- 2 T 1 1 2 1 4 17 6 

tan - = 1 + — x z + — T 4 + x b + 



(18) 



x 2 12 120 20160 

the exact evolution is reproduced, provided the timestep is restricted to the interval 
< x < it. If we interchange the roles of //kick and // m0 ve one combination instead 
gives 

/s i i 1 — jinkx 2 (1 — ^mk)mx 



-kx 1 — \kmx 2 




(19) 
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which becomes exact if we choose 



2 x , suit 
m = - tan-, k= , (20) 

T 2 X 

again provided the timestep is restricted to the interval < x < n. 

2.2. Higher-dimensional linear systems 

It should be clear that this idea works for systems of harmonic oscillators in 
general, i.e. for quadratic Hamiltonians of the form 

H(q, P )= l -{p T M P + q T Kq), (21) 

where M and K are symmetric matrices. For a choosen splitting scheme and step 
interval x there are always modified matrices, 

X 2 X 4 X 6 

M T =M — —M (KM) + —M (KM) 2 - — M (KM) 3 + ^(t 8 ), (22) 
K % = K+^(KM)K+^(KM) 2 K+^^(KM) 3 K+0(x s ) 1 (23) 

generating kick-move-kick flow which reproduces the exact flow up to corrections 
of order T 8 . It should be obvious how this can be extended to arbitrary order in x , 
with coefficients taken from the expansions in equation ( [T8] ). In principle this can 
be used to reproduce the exact flow, provided x is not too large. The constraint is 
the restriction < x < nj tt) m ax> where (Omax is the largest frequency of the system. 

3. General potentials 

For a more general treatment we consider Hamiltonians of the form 

H(q, P )= l -p T M P + V(q). (24) 
A series solution of the Hamilton equations in powers of x is 

q a e =q a + p a X - \d a Vx 2 - ld a (DV)x 3 + 0(x 4 ), 
2 6 

Pi = Pa ~ d a Vx - l -d a (DV)x 2 + d a f^DV - l -D 2 v\ x 3 + &(x 4 ). 
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Here we have introduced notation to shorten expressions, 

da^^~ a , d a = M ab d b , p a =M ab p b , D = p a d a , D=(d a V)d a , 

where we employ the Einstein summation convention: An index which occur 
twice, once in lower position and once in upper position, are implicitly summed 
over all available values. I.e, M ah d\, = Y,bM ab di ? (we will generally use the matrix 
M to rise an index from lower to upper position). 

If we instead use a splitting method to generate the flow, with generators H\ = 
V(q) and H 2 = \p J Mp = T (i.e., a kick-move -kick scheme), we obtain 

q a = q a + p a X - X -d a V X 2 + ^(T 4 ), 

Pa = Pa ~ d a VT ~ \d a {DV)% 2 + d a Qw - ^D 2 V^j T 3 + 0{% A ). (25) 

As expected the result differs from the exact result in the third order. However, 
the difference can be corrected by modifiying the generators, H\ — > T + T% and 
H 2 -t-V + V 2 , with 

T 2 = -^D 2 Vt 2 , V 2 = ^DVx 2 . (26) 

Specialized to a one-dimensional system with potential V = \q 2 this agrees with 
equation ( fl"8| ). With this correction the kick-move-kick splitting scheme agrees 
with the exact solution to 4 th order in T, but differ in the T 5 -terms. We may again 
correct the difference by introducing fourth order generators, H\ — >■ T + T 2 + T4 
and H 2 -> V + V 2 + V 4 , with 

74 = 720 ^ ~ 9D ° 2 + 3DDD ^ V Va = (27) 
Specialized to a one-dimensional system with potential V = \q 2 this agrees with 
equation ( [T8] ). With this correction the kick-move-kick splitting scheme agrees 
with the exact solution to 6 th order in T, but differ in the T 7 -terms. We finally 
correct this difference by introducing sixth order generators, U\ — > T + T 2 + T4 + 
T 6 and H 2 -> V + V 2 + V 4 + V 6 , with 

T 6 = - (2D 6 - AODD 4 + A6DDD 3 - \5D 2 DD 2 + 5AD 2 D 2 

- 9 DDDD - 42 DD 2 D + 1 2 D 2 D 2 ) V T 6 

(28) 
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where we have introduced 



D 3 = (d a V)(d b V)(d c V)d a d b d c . (29) 

Specialized to a one-dimensional system with potential V = \q 2 this agrees with 
equation < [T8"j ). With this correction the kick-move-kick splitting scheme agrees 
with the exact solution to 8 th order in z, but differ in the T 9 -terms. The process 
may be continued to higher orders in z, 

Hl ^ T+ ^T 2k , H 2 ^V+Y,V 2k . (30) 

k>\ k>\ 

To keep track of the algebraic expressions which occured during the calcu- 
lations above, we have represented them graphically in terms of bi-colored tree- 
diagrams. I.e., these calculations are related to "rooted-tree-type" theories. Our 
tree-diagrams describing T 2 k and V 2 k, and the generating functions below, are 
unrooted (the derivatives of these scalar functions can be represented by rooted 
trees). It is fairly straightforward to find the general structure of the order z N cor- 
rection terms, but more laborious to compute the rational coefficients multiply- 
ing each term. They are simplest found by considering enough special cases for 
unique determination. After the explicit expressions (26} 27]) were found we ver- 



ified them manually for a general Hamiltonian (21 ) using graphical calculations. 



The explicit expression (28 ) has been checked against a general Hamiltonian (21 ) 



acting on a four-dimensional phase space (i.e., with two-dimensional q and p). 



4. Solving the move steps 

Addition of extra potential terms V — > V e ff = V + V 2 + V4 + . . . , is in principle 
unproblematic for solution of the kick steps. The equations, 

<f = 0, p a = -d a V eS (q), (31) 

can still be integrated exactly, preserving the symplectic structure. The situation 

is different for the kinectic term T — > T e ff = T + T 2 + T4~\ , since it now leads 

to equations 

q a = 3— T eS (q,p), Pa = -d a T eS (q,p), (32) 
a Pa 

which are no longer straightforward to integrate exactly. Although the problem- 
atic terms are small one should make sure that the move steps preserve the sym- 
plectic structure. Let q,p denote the positions and momenta just before the move 
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step, and Q,P the positions and momenta just after. The relation between q,p 
and Q,P can be expressed in terms of a generating function (cf. Golstein [|22|. 
Arnold ll23l), 

G(g, P; t) = q a P a + AG(q, P; x) , (33) 
such that the transformation 

ff = w» »" = W" =P " + W <34) 

preserves the symplectic structure exactly. However, note that the relation be- 
tween p and P in general is a nonlinear equation of the form 

P a = Pa-j ll AG(q,P;x), (35) 

where the second term on the right is of order T 3 or higher. We solve this equation 
by iteration. With p(°) = p, 



Pa -—AG(q,P^;x). (36) 



Writing = P + APW, with P the exact solution, we find to first order in 
AP that 

Ap!/ i+1) = -j^AG{q,P;T)Alf> = -AG fl *AP« (37) 

Let A be the eigenvalue of AG a b with largest magnitude. Then the iteration con- 
verges exponentially fast towards the exact solution, with AP^ decaying like 
X n ~ t 3 ' 1 . Since it is most to gain by a higher order method when the timestep x 
is small, we assume that A is small in cases of practical relevance. Our experi- 
ence is that the iteration scheme is robust, with 3-4 iterations been sufficient for 
computations to double precision accuracy. 

Some of our theoretical results have already been given in the literature. The 
generating function formalism has been used earlier by Feng [fT8l and Feng et.al 
[fT9l to construct canonical difference schemes (see also Channell and Scovel ||20l , 
Stuchi BTlO . They give the result ([26]), but the actual solution of the resulting 



implicit equations are not discussed. One can construct a generating function 
for the full symplectic evolution over a timestep x, without combination with a 
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splitting method. However, in that case the resulting nonlinear equations would 
be more time consuming and/or difficult to solve by direct iteration. 

We now explicitly construct G so that the move step is reproduced to sufficient 
accuracy. Consider first the case when H\ = T . The choice G = q a P a + \P a P a x 
gives 

Q a = q a +P a T, Pa=Pa, 08) 

which is the correct relation. Now add the ^-term to the move step. The exact 
solution of equation (32) becomes 

Q a = q a + p a T- \d a DV X 3 - ±d a D 2 V T 4 + 0(x 5 ), 



Pa = Pa + T2d a D 2 V X 3 + ±d a D 3 V + 6{% 



(39) 



Compare this with the result of changing 



G^G- 



12" 



(40) 



where @ = P a d a . The solution of equation (34) change from the relations (38 ) to 



Pa 



q a + pa % _ l d a@y T 3 _ l^^ly % * + ff ^5 



Pa 



12' 



l 2 Vz 3 -±d a @ 3 VT 4 + 



(41) 
(42) 



Since Qi is linear in P, equation (42) constitute a system of third order algebraic 
equations which in general must be solved numerically. This should usually be a 
fast process for small x. An exact solution of this equation is required to preserve 
the symplectic structure, but this solution should also agree with the exact solution 
of (32) to order T 4 . This may be verified by perturbation expansion in T. A 
perturbative solution of equation (42) is 

Pa = Pa + Ty^VX 3 + ld a D 3 VZ 4 + ^(T 5 ), 



which inserted into (41 ) reproduces the full solution (39) to order t . 

This process can be systematically continued to higher orders. We write the 
transformation function as 

CO 



G T 



EG*** 

k=0 



(43) 
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and find the first terms in the expansion to be 

G = q a P a , G x = \P a P Ul G 2 = 0, G 3 = -±@ 2 V, G 4 = -^ 3 V, 
G 5 = -2^5 (3@ 4 + 3D@ 2 -%iD@) V, 
G 6 = -^ K) {2 2> 5 +8D@ 3 - 5 ^D5> 2 ) V, 

G 7 = -20T6O (lO# + 10D# + 90^D#-75^ 2 D^ 2 

+ 18D 2 ^ 2 -3D^D^-14^D 2 ^ + 45> 2 D 2 )V, 
Gg = - 45^20 (3 2) 1 - 87 D3> 5 + 23 1 9D9 A - 1 33 2i 2 D3> 3 + 63 D 2 ^ 3 - 3 ^D 2 ^ 2 - 2 1 ^ 2 D 2 9 

+ 4 #D 2 - 63 + 25 V. 

Also in these calculations we represent the algebraic expressions by bi-colored 
tree diagrams, to better visualise and understand their structure. The possible 
graphical structures for G n is fairly simple to write down. But it is quite laborious 
to find the rational coefficients multiplying each graph. They are simplest found 
by considering enough special cases for unique determination. After that we have 
verified the expressions up to manually using graphical calculations, and G7, 
G% against a general Hamiltonian ( pi) acting on a four-dimensional phase space 
(i.e., with two-dimensional q and p) using a computer algebra program. 

5. Numerical results on nonlinear systems 

5.1. One-dimensional anharmonic oscillator 

It remains to demonstrate that our algorithms can be applied to real examples. 
We have considered the Hamiltonian 

H= 1 -p 2 + 1 -q\ (44) 

with initial condition q(0) = 0, p(0) = 1. The exact motion is a nonlinear oscilla- 
tion with H constant equal to \, and period 

T = 4 f 2 J^L = 2 1 /4 B( 1 ,1)^6.236339.... (45) 
•A) yj2-q A 4 2 

Here B(;c,y) = T(x)T{y) /T(x + y) is the beta function. In Figure [T] we plot the 
behaviour of (H —I) / z 2+n during the last half of the 16 th oscillation, for various 
values of T and corrected generators up to order T 6 (corresponding to n = 6). 
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Figure 1: This figure illustrate how well energy is conserved with the various splitting schemes. 
The quanties plotted is (H — ^ ) / z m for T = 0.2 (squares), T = 0. 1 (triangles) and T = 0.05 (lines). 
Here m — 2 for the Stormer-Verlet scheme (dotted line), m = 4 for the T 2 -corrected generators 
(dash-dotted line), m = 6 for the T 4 -corrected generators (dashed line), and m = 8 for the T 6 - 
corrected generators (fulldrawn line). Each plotted quantity is essentially the value of the next 
correction at the visited point in phase space. Since the plot is taken over the last half of the 16 th 
period the figure also give some indication of how well the exact oscillation period is reproduced 
by the scheme. The deviation is quite large for the Stormer-Verlet scheme when T = 0.2; to avoid 
cluttering the figure we have not included these points. 



5.2. Fermi-Pasta-Ulam-Tsingou type problems 

Here we will consider a one-dimensional closed chain of d particles (as illus- 
trated in Figure [3]) interacting with its nearest neighbours through a potential U, 
and possibly with a local substrate through a potential V. The latter will confine 
the n th particle to the vicinity of a position = nL/d, where L is the circumfer- 
ence of the chain. 

The class of models for this system include the Fermi-Pasta-Ulam-Tsingou 
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Energy preservation with the Stormer-Verlet splitting scheme Energy preservation with ^-'-corrected generators 




25595 25596 25597 25598 25599 25000 1038395 1638396 1038397 1038398 1038399 1638400 



Figure 2: These figures illustrate the long time behaviour through the last half of the 16 period 
for the Stormer-Verlet, first half of 257 th period for T 2 -corrected, last half of the 4 104 th period for 
T 4 -corrected and last half of the 262718 th period for T 6 -corrected schemes. Different timesteps T 
have an effect on the period of oscillation, but the preservation of energy remains stable for a very 
long time. 

(FPU) problem introduced in 1953 by Fermi et. al. Il24l for investigating equipar- 
tition of energy among the different degrees of freedom. This study opened up 
the doors of research in many fields of mathematics and physics. Lot of research 
is going on in different field of studies to understand the highly unexpected dy- 
namical behaviors. A preview of last 50 years comprehensive study on FPU is 
provided by G.P Berman [|25l . 

In the recent papers by E. Hairer and C. Lubich ll26l presented a study for 
FPU by using the modulated Fourier expansions on the chains with large number 
of particles, and in ll27l I. McLachlan and D. Neal have given the good comparison 
of integrators on FPU. A good analysis of the Stormer-Verlet method applied to 
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the FPU problem has been given by G. Benettin and A. Ponno [|28ll by using BCH 
formula. Numerical methods with some results on FPU are the part of ll29l . 

Here we will demonstrate that our 
integrators can be implemented and 
applied in practise to these type of 
models. There is of course a compu- 
tational cost per timestep by going to 
a higher order method, but asymptoti- 
cally that cost grows linearely with the 
size of the system. There is also a cost 
in complexity of code implementation, 
which we have solved by writing a pro- Figure 3 . Chain of neighbouring particles 
gram for automatic generation of the 
numerical code l30l . 

Let q m {t) = r m {t) —t% m , where r m (t) the position of the n th particle, and con- 
sider the system described by the Hamiltonian 

# (q,P) = I £ Pl + £ V{9m) + (46) 
m—Q m=0 m=0 

where d is the number of particles, and s m = q,„+i — q m . A class of model which 
includes both the linear chain and the FPU model can be obtained by choosing 

V(q)= l -OJ 2 q\ U(s)= l -s 2 + ^as 3 + ^s 4 , 

where U is describing the interactions between particles. This model has been 
referred to as the FPU a + /3 model (with CO 2 = 0). In this example we have 
used a = and j8 = 1. We have tested the methods with respect to (i) energy 
conservation, (ii) deviation of the generated solution from the exact solution^ and 
(iii) efficiency of the methods with respect to CPU time. 

In Figure|4]we show the scaled energy error on FPU by using different choices 
of T of all four methods. For these experiments we consider 9 particles with 
initial energy £'(0) = 1.425. As can be seen the energy conserved very well for all 
methods, with the error scaling like for a method of order N. As demonstrated 
by the long time behaviour in Figure[5]the energy error does not increase noticably 
with time. 



Actually a numerical solution of the same system computed to very high precision. 
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Figure 4: Scaled energy error for higher order methods of the FPU 



Another quantity of interest in a system with many degrees of freedom is the 
global error, i.e. a measure how much the numerical solution deviates from the 
exact solution. Here an exact solution is not available. Instead we have generated 
a very accurate solution by use of our 8 th order method with timestep x = 5 ■ 10 -4 , 
calculated with multiprecision (50 decimal digits) floating point accuracy. This is 
for practical purposes as good as an exact result, and we will refer to it as such. 

We have investigated several measures of deviation; they all give qualitatively 
the same results. Here we will only discuss the quantity 



£(t) = \\(q(t),p(t))-(Qn,Pn)\\ 2 



d-l 

E 

m—0 



y, (flmif) Gm,n) Cp#»(0 Pm.n 



1/2 



(47) 
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Long time scaled energy error 



2nd order 



4th order 
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Figure 5: Long time scaled energy error for a Fermi-Pasta-Ulam-Tsingou type problem, computed 
with the Stormer-Verlet (A) and higher order corrected integrators (B-D) with timestep T = ^ 
and initial energy E(0) = 1.425. 



where q n (p n ) denote the positions (momenta) of the numerical solution at a 
timestep n such that nx = t, and q(t) ip{t)) denote the positions (momenta) of the 
exact solution at time t. As shown in Figure [6] the global error behaves roughly 
like 

£{t)~Ctx N , (48) 

for relatively short times t. Here C is a constant which depends on the order iV of 
the method and the initial conditions. This is in agreement with exact behaviour 
of integrable systems, cf. Theorem 3.1 in the book [Q]| by Hairer et. al. 
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Global error at a fixed time Global error versus time 




1(T 2 KT 1 5 10 15 20 



Figure 6: The left frame shows how the global error (here measured at time t = 10) depends on the 
timestep T and the order N of the integration scheme. As expected this error varies like T N (as long 
as it is small). The right frame shows how the global error grows with time, here for a timestep T = 
4^. The dashed lines are eyeball fits to linear error growth, cf. e(f) = || (q(t),p(t)) — (<7„,p„)|| 2 ~ 
Ct Tr . These results are for a lattice of d = 9 particles. 



To check the efficiency of our methods in practical use, we have also measured 
CPU time used to integrate systems with different number d of particles, with 
d ranging from 9 to 50000. All runs have been done on the same system, a 
workstation equipped with two six-core Opteron 2431 processors, but using code 
written in NumPy. Hence, the code is not parallellized and run on a single core. 
Some results, run with timestep % = 1/12 for all methods, is shown in Figure [7] 
Under these conditions we find that the CPU time increases by a factor of about 
10 for each step in order. From the left frame of Figure [6] we see that this step also 
increases the accuracy with a factor of about 10 _1 T 2 (for d = 9 particles). If we 
want a prescribed accuracy I0~ p for the global error e(t) at time t we may choose 
to use lower order method with a small timestep (which requires many steps n), or 
a higher order method with fewer, but more time-consuming steps. Which choice 
is best? For the parameters displayed in Figure [6] we estimate the condition 

e(t)*10- 2 - N / 2 tx N *10- p . (49) 

I.e., we must choose a timestep such that 

x N « - x lO 2 ^/ 2 ^, (50) 
t 

which requires 

n ~ I ~ _L x t (l+l/A0 io(p-V/n (51) 
T v^To 
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steps, where each steps requires a CPU time ? step ~ to 10^ for some constant to 
which depends on the computer being used. Hence, we should choose N to mini- 
mize 

Thro = "'step « -^L x lO^^^io*-^. (52 ) 



10 

Treating N as a continuous varible gives the optimal value 



iV p t ~ V 2 ( p + lo gio ? -2)- 



(53) 
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Figure 7: CPU time 7cpu used to solve a lattice of d particles for 1000 timesteps (t = -pj) for 
schemes of different orders N. Asymptotically, 7cpu grows linearly with d. The penalty for 
increasing the order N by 2 is about a factor 10 increase in 7cpu, when d and the number of 
timesteps is kept fixed. 



6. Concluding remarks 

In this paper we have shown that it is possible to systematically extend the 
standard Stormer-Verlet symplectic integration scheme to higher orders of accu- 
racy, and that the higher order schemes can be applied in practise to physical 
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systems of interest, including FPU-like lattice problems with many particles (with 
nearest-neighbour interactions). As illustrated by equation p3| ), it is advantageous 
to use a higher order method when one wants a solution of high precision P, and 
also if one wants a solution of moderate accuracy but over a long time interval. 

As demonstrated, the theoretical algorithms have been implemented and tested. 
One rapidly discovers that it is a nightmare to do a correct implementation by 
hand. The general compact form of these schemes usully expand to very long 
expressions, which are laborious and error-prone to handle manually. We have 
therefore developed a set of computer routines which automatically generate the 
basic numerical integrators for a complete timestep of each specific model. 

For the cases we have investigated these integrators perform according to ex- 
pectations, sometimes even better than expected. 
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